{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# B 13-2 Logistic 回归的参数分析和假设检验\n",
    "\n",
    "结果变量为二项分类变量，此时不适用线性回归，应使用 Logistic 回归。\n",
    "\n",
    "## 案例\n",
    "\n",
    "研究吸烟 (X1)、饮酒 (X2) 与食管癌 (Y) 之间的关系。\n",
    "\n",
    "### 数据"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>X1</th>\n",
       "      <th>X2</th>\n",
       "      <th>Y</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>881</th>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>882</th>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>883</th>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>884</th>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>885</th>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>886 rows × 3 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "     X1  X2  Y\n",
       "0     0   0  1\n",
       "1     0   0  1\n",
       "2     0   0  1\n",
       "3     0   0  1\n",
       "4     0   0  1\n",
       "..   ..  .. ..\n",
       "881   1   1  0\n",
       "882   1   1  0\n",
       "883   1   1  0\n",
       "884   1   1  0\n",
       "885   1   1  0\n",
       "\n",
       "[886 rows x 3 columns]"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import polars as pl\n",
    "\n",
    "with open(\"B_13_2-data.csv\", \"r\") as reader:\n",
    "    header = reader.readline().split(\",\")\n",
    "\n",
    "    df = pl.DataFrame(schema={h: pl.Int8 for h in header[:-1]})\n",
    "\n",
    "    while line := reader.readline():\n",
    "        line = line.strip().split(\",\")\n",
    "\n",
    "        for i in range(int(line[-1])):\n",
    "            new_df = pl.DataFrame(data=[line[:-1]], schema={h: pl.Int8 for h in header[:-1]}, orient=\"row\")\n",
    "            df = df.vstack(new_df)\n",
    "\n",
    "df = df.to_pandas(use_pyarrow_extension_array=False)\n",
    "\n",
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Logistic 回归"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization terminated successfully.\n",
      "         Current function value: 0.654301\n",
      "         Iterations 5\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>Logit Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>           <td>Y</td>        <th>  No. Observations:  </th>  <td>   886</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                 <td>Logit</td>      <th>  Df Residuals:      </th>  <td>   883</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>                 <td>MLE</td>       <th>  Df Model:          </th>  <td>     2</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>            <td>Mon, 09 Dec 2024</td> <th>  Pseudo R-squ.:     </th>  <td>0.05582</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                <td>22:04:36</td>     <th>  Log-Likelihood:    </th> <td> -579.71</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>converged:</th>             <td>True</td>       <th>  LL-Null:           </th> <td> -613.98</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>     <td>nonrobust</td>    <th>  LLR p-value:       </th> <td>1.305e-15</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "       <td></td>         <th>coef</th>     <th>std err</th>      <th>z</th>      <th>P>|z|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Intercept</th>  <td>   -0.9099</td> <td>    0.136</td> <td>   -6.699</td> <td> 0.000</td> <td>   -1.176</td> <td>   -0.644</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(X1)[T.1]</th> <td>    0.8856</td> <td>    0.150</td> <td>    5.904</td> <td> 0.000</td> <td>    0.592</td> <td>    1.180</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(X2)[T.1]</th> <td>    0.5261</td> <td>    0.157</td> <td>    3.348</td> <td> 0.001</td> <td>    0.218</td> <td>    0.834</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/latex": [
       "\\begin{center}\n",
       "\\begin{tabular}{lclc}\n",
       "\\toprule\n",
       "\\textbf{Dep. Variable:}   &        Y         & \\textbf{  No. Observations:  } &      886    \\\\\n",
       "\\textbf{Model:}           &      Logit       & \\textbf{  Df Residuals:      } &      883    \\\\\n",
       "\\textbf{Method:}          &       MLE        & \\textbf{  Df Model:          } &        2    \\\\\n",
       "\\textbf{Date:}            & Mon, 09 Dec 2024 & \\textbf{  Pseudo R-squ.:     } &  0.05582    \\\\\n",
       "\\textbf{Time:}            &     22:04:36     & \\textbf{  Log-Likelihood:    } &   -579.71   \\\\\n",
       "\\textbf{converged:}       &       True       & \\textbf{  LL-Null:           } &   -613.98   \\\\\n",
       "\\textbf{Covariance Type:} &    nonrobust     & \\textbf{  LLR p-value:       } & 1.305e-15   \\\\\n",
       "\\bottomrule\n",
       "\\end{tabular}\n",
       "\\begin{tabular}{lcccccc}\n",
       "                    & \\textbf{coef} & \\textbf{std err} & \\textbf{z} & \\textbf{P$> |$z$|$} & \\textbf{[0.025} & \\textbf{0.975]}  \\\\\n",
       "\\midrule\n",
       "\\textbf{Intercept}  &      -0.9099  &        0.136     &    -6.699  &         0.000        &       -1.176    &       -0.644     \\\\\n",
       "\\textbf{C(X1)[T.1]} &       0.8856  &        0.150     &     5.904  &         0.000        &        0.592    &        1.180     \\\\\n",
       "\\textbf{C(X2)[T.1]} &       0.5261  &        0.157     &     3.348  &         0.001        &        0.218    &        0.834     \\\\\n",
       "\\bottomrule\n",
       "\\end{tabular}\n",
       "%\\caption{Logit Regression Results}\n",
       "\\end{center}"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                           Logit Regression Results                           \n",
       "==============================================================================\n",
       "Dep. Variable:                      Y   No. Observations:                  886\n",
       "Model:                          Logit   Df Residuals:                      883\n",
       "Method:                           MLE   Df Model:                            2\n",
       "Date:                Mon, 09 Dec 2024   Pseudo R-squ.:                 0.05582\n",
       "Time:                        22:04:36   Log-Likelihood:                -579.71\n",
       "converged:                       True   LL-Null:                       -613.98\n",
       "Covariance Type:            nonrobust   LLR p-value:                 1.305e-15\n",
       "==============================================================================\n",
       "                 coef    std err          z      P>|z|      [0.025      0.975]\n",
       "------------------------------------------------------------------------------\n",
       "Intercept     -0.9099      0.136     -6.699      0.000      -1.176      -0.644\n",
       "C(X1)[T.1]     0.8856      0.150      5.904      0.000       0.592       1.180\n",
       "C(X2)[T.1]     0.5261      0.157      3.348      0.001       0.218       0.834\n",
       "==============================================================================\n",
       "\"\"\""
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import statsmodels.formula.api as smf\n",
    "\n",
    "# perform Logistic regression\n",
    "\n",
    "# define the formula (response ~~ predictors)\n",
    "formula = \"Y ~ C(X1) + C(X2)\"\n",
    "\n",
    "# fit the model\n",
    "model = smf.logit(formula=formula, data=df).fit()\n",
    "\n",
    "model.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Logistic 回归结果显示，常数项的系数为 0.136，吸烟 (X1)、饮酒 (X2) 的回归系数分别为 0.150 和 0.157。 $ \\forall \\ p < 0.05 $ ，说明喝酒和吸烟对回归方程都有统计学意义。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在此次 Logistic 回归中，X1, X2 都是二项分类变量，不存在单位上的差异。不知道课本上所说的“标准化回归系数”有何意义。"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": ".venv",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
